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Over the past few years Fe chalcogenides (FeSe/Te) have advanced to the forefront of Fe-based 
superconductors (FeBS) research. The most intriguing results thus far are for intercalated and monolayer 
FeSe, however experimental studies are still inconclusive. Yet, bulk FeSe itself remains an unusual case 
when compared with pnictogen-based FeBS, and may hold clues to understanding the more exotic FeSe- 
derivatives. The FeSe phase diagram is unlike the pnictides: the orthorhombic distortion, which is likely 
to be of a "spin-nematic" nature in numerous pnictides, is not accompanied by magnetic order in FeSe, 
and the superconducting transition temperature T c rises significantly with pressure before decreasing. In 
this paper we show that the magnetic interactions in chalcogenides, as opposed to pnictides, demonstrate 
unusual (and unanticipated) frustration, which suppresses magnetic, but not nematic order, favors ferro- 
orbital order in the nematic phase and can naturally explain the nonmonotonic pressure dependence of 


the superconducting critical temperature T C (P). 

While full consensus regarding the mechanism of high- 
temperature superconductivity in Fe-based superconductors 
(FeBS) remains elusive, nearly all researchers agree that it 
is unconventional and that it has a magnetic origin 1 2 . How¬ 
ever, there is a divergence of opinions on the nature of the 
electrons responsible for magnetism. There is an itinerant 
approach based on calculating the spin susceptibility with 
moderate Coulomb (Hubbard) and Hund’s interactions 3 10 
as well as a localized approach where itinerant electrons re¬ 
sponsible for conduction and the Fermi surface interact with 
local spin^MEH. Finally, there is an increasingly popular de¬ 
scription where the electrons have a dual character and pro¬ 
vide the local moments, the interaction between them, and 
the electronic conductivity 1 - 3 '^. Within this picture, FeBS 
can still be reasonably mapped onto a short-range model of 
pairwise interactions between the local moments. 

Following the discovery of the FeBS, there were multi¬ 
ple attempts to map the exchange interactions onto the 
Heisenberg model. The J 1 -J 2 model on the square lat¬ 
tice 17 with nearest- (Ji) and next-nearest-neighbor (J 2 ) ex¬ 
change couplings was a natural starting point®*®, but re¬ 
quired dramatically different couplings for ferro- and antifer¬ 
romagnetic neighbors, J\ a <C Jib to reproduce the observed 
spin wave^®®and ab-initio calculations!® it also failed to 
describe the double-stripe configuration in FeT d— I— 4 The 
model was extended to include third-neighbor exchange J 3 27 
to reproduce the FeTe magnetic ground state. However, 
only the Ising model has this configuration as a solution, 
and in the Heisenberg model it is not a ground state for any 
set of parameters 28 29 . Therefore adding J 3 does not solve 
the problem. Besides, the Ji a <C Jib implies an unphysical 
temperature dependence of the exchange constants (as T 
approaches Tn, by symmetry J\ a —1 Jib). 

There were attempts to overcome these problems by 
adding the nearest-neighbor biquadratic exchange interac¬ 
tion K(Si • Sj ) 2 to the Ji-J 2-^ 3 ®- 31 or J 1 -J 2 -J 3 32 Heisen¬ 
berg model. The three-neighbor Heisenberg model with bi¬ 


quadratic term (denoted J 1 -J 2 -J 3 -K model from now on) 
eliminates the need for the J\ a ,ib anisotropy of the nearest- 
neighbor exchange and, for sufficiently large K and J 3 , has a 
ground state consistent with that of FeTe. The biquadratic 
coupling in this model is also essential to explain the split¬ 
ting between the antiferromagnetic and orthorhombic phase 
transitions in the Fe pnictideJ® 33!! ®. 

Whereas the magnetism in Fe-pnictides is successfully ex¬ 
plained by the J 1 -J 2 -J 3 -A' model, the Fe-chalcogenides re¬ 
main problematic. Specifically, there are two important un¬ 
resolved controversies regarding bulk FeSe; (i) it shows a 
structural transition at T s ~ 90 K but, contrary to the Fe- 
pnictides, no magnetic order is observed below T s . Instead, 
an extended nematic region is detectec 35 -*® and the system 
becomes superconducting at T c ~ 8 K. (ii) The supercon¬ 
ducting T c first increases with pressure and then decreases, 
forming a dome 37 . This is in apparent contradiction with 
the expectation of a decreasing T c with pressure when mag¬ 
netism is absent. 

In the present work we propose a solution to this mystery 
and generalize the results to the family of Fe-chalcogenides 
FeSe/Te. We show, using ab-initio density functional theory 
calculations and effective model considerations, that many 
properties of FeSe/Te are related to its unusual magnetic 
frustration, absent in the Fe-pnictides. We show that J1-J2- 
J 3 -K is the minimal spin model that includes the relevant 
complexity of the magnetism in Fe-chalcogenides. We then 
identify a new range of parameters appropriate for FeSe/Te 
where a highly competitive novel “staggered dimer” phasd 3 ^ 
is stabilized (recently shown to be the ground state of FeSe 
in ab-initio calculation^). 


Exchange model and phase diagram 
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We define the J1-J2-J3-K model on the square lattice as 


H = 22 ■ liij + K (m, • m ^) 2 
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(i) 


The first sum is taken over all nearest neighbor {i,j} pairs 
of Fe spins, the second one over all next nearest neighbors, 
etc. m is the unit vector in the spin direction, |m,| = 1. 
For K = 0(oo), this model reduces to the already solved 
Heisenberg® (Isinf!®) on the square lattice. 

To begin, we review the phase diagrams for the standard 
J1-J2-J3 Ising and Heisenberg models. In Fig. [2^ we show 
the mean-field T = 0 Ising phase diagram which includes 
the staggered dimer and double stripe ground states (see 
Fig. [T| for pattern definition). While this phase diagram 
can explain the ab-initio magnetic states of FeSe and FeTe, 
note that the Ising model is inapplicable to low-anisotropy 
materials such as the pnictides and chalcogenides, and the 
Heisenberg model is more appropriat e! 24 * 30 ^ 2 . The mean- 
field phase diagram for the Heisenberg model at T = 0 
is shown in Fig. [ 2 J 3 (quantum corrections introduce minor 
changed®). The double stripe phase has measure zero [it 
is a degenerate case of q = ( Q,Q )]■ This means that no 
Heisenberg model can explain the formation of a colli near 
double stripe state. 

The review of the Ising and Heisenberg phase diagrams 
elucidates the two theoretical problems that have been un¬ 
deremphasized in previous analyses of the magnetic inter¬ 
actions of the Fe-based superconductors, especially in the 
chalcogenides: (1) the Heisenberg model does not account 
for all relevant magnetically ordered states and, by impli¬ 
cation, does not properly describe spin fluctuations, and 
(2) the single and double stripe magnetic states are not 
the only important ground state candidates for the chalco¬ 
genides; there is a third one, the staggered dimers, which 
is highly competitive, but has been routinely ignored. To 
address these problems the biquadratic term, K, needs to 
be quantitatively taken into account. 

We solved the full J 1 -J 2 -J 3 -K model (Equation (JT])) for 
general K in the mean-field limit and found six possible 
ground states (see Supplementary Materials). Hu et. a/!® 
attempted prevously to solve this model, but missed the 
staggered dimer phase® which, we argue, is the key to un¬ 
derstanding FeSe. A representative example phase diagram 
with K = 0.1 is shown in Fig. [2J:. For a small, but non-zero 
K and J3 the staggered dimer phase becomes stable in a 
narrow (| J 2 — Ji/2| < interval near the critical 

value J\ = 2 J 2 , and at sufficiently large J 3 (J 3 > J 2 /8A') 
the collinear double stripe structure is stabilized. As I\ 
grows, these collinear regions also grow, and at K > Ji/2 
the phase diagram becomes identical to the Ising phase di¬ 
agram in Fig. [2^. Actual materials will be seen to lie in the 
intermediate region, 0 < K < Ji/2. 

First-principles calculations 

We performed density functional theory (DFT) calcula¬ 
tions to obtain parameters for Equation ([I]) and place FeSe 


and FeTe into the context of the J1-J2-J3-K phase dia¬ 
gram. There is a caveat though: due to the itinerant char¬ 
acter of magnetism in FeBS, mapping onto local moments 
models such as Equation (jT|) has limited accuracy. A fun¬ 
damental assumption of the standard Heisenberg model is 
that the magnetic moments are rigid, and this is an excel¬ 
lent assumption for systems with highly localized electrons, 
such as the high-T c cuprates, but relatively poor for itiner¬ 
ant electrons. Magnetic interactions in metals tend to have 
long range tails, non-pairwise interactions, and the moments 
may depend on the magnetic ordering pattern. A clear ex¬ 
ample of the failure of the Heisenberg-biquadratic models is 
that the double stripe (Fig. [TJd) and plaquette (see Fig. Ip 
in the Supplementary Material) configurations are degener¬ 
ate in any such model, but in DFT the double stripe is 8 
meV/Fe lower in energy than the plaquette configuration®. 
Therefore, we cannot expect to derive a Heisenberg model, 
with or without the biquadratic K, that captures exactly 
the energetics of all possible magnetic configurations. 

Despite these limitations, the J3-J2-J3-K model is the 
simplest framework that accounts for all the magnetic 
ground states that DFT and experiment find in different 
FeBS, and arguably is also the most complex one that still 
allows for an analytic solution. Since we are interested in 
spin fluctuation-driven effects such as superconductivity and 
spin-nematicity, which are low-energy phenomena, we es¬ 
tablish a set of criteria for our fits, with the main goal to 
select a consistent set of magnetic states and obtain pa¬ 
rameters that reproduce the low-energy hierarchy obtained 
within DFT. The criteria are detailed in the Methods sec¬ 
tion, and the chosen magnetic structures are shown in Fig. [l] 
panels a through e. 

We performed calculations for FeSe at three representa¬ 
tive pressures of 0, 4, and 9 GPa, and for FeTe at ambient 
pressure, see the Methods and Supplementary Materials for 
details. In all cases we used experimental lattice and in¬ 
ternal parameters in tetragonal structures, as discussed in 
Methods. We fitted to the five magnetic configurations 
reported in Fig. [3] and extracted the J\, J2, and J 3 parame¬ 
ters. The biquadratic term was extracted from noncollinear 
calculations as in Ref.351 The resulting J 1 -J 2 -J 3 -K model 
parameters are reported in Table [Tj Note that the error bars 
reflect the fit inaccuracy, and not the much smaller errors 
of the underlying DFT calculations. 

First of all, we confirmed that the "staggered dimer” 
configuration!® is 13 meV/Fe lower in energy than the sin¬ 
gle stripe configuration and is the true DFT ground state 
for FeSd 39 (see Fig. [3]). The same phase is also the lowest 
in energy in FeTe, as long as one does not take into ac¬ 
count the magnetoelastic coupling. The calculated energy 
difference between the double stripe and staggered dimer 
configurations in tetragonal FeTe is tiny, ~ 1 — 2 meV/Fe. 
However, upon full structural relaxation into a monoclinic 
structure the double stripe pattern gains more magnetoe¬ 
lastic energy than the staggered dimer one (which relaxes 
into an orthorhombic structure) and ends up lower by a few 
meV, with the crystallographic distortion in agreement with 
experiment 25 ^®. 
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Another important result is that while the main con¬ 
tenders for the ground state of FeTe are the double stripe 
(fids = (7 t/2, 7t/2)) and the staggered dimer (q^ = 
(7r, 7t/2)) structures, with the staggered trimers (q tr i = 
( 7 r, 7 r/ 3 )) a close third, in FeSe the double stripe structure 
is not competitive at all. In FeSe the lowest energy states 
are the staggered dimers, trimers, tetramers and single 
stripes, with q, respectively, (7t,7t/2), (n, 7t/3), (7t,7t/4), 
and (q ss = (n, 0)). From this, one can conclude that while 
in experiment the long range order of FeSe is destroyed by 
spin fluctuations, the most relevant ones are those with the 
corresponding wave vectors as listed above, and, very likely, 
with any q = (- 7 T, Q) such that 0 < Q < 7 r/ 2 . 

Importantly, when FeSe is structurally optimized in any of 
the low energy magnetic structures, it admits an orthorhom¬ 
bic structure quantitatively consistent with the experiment, 
(a — b)/(a + b) ~ 0 . 2 %, while optimization without mag¬ 
netism never breaks the tetragonal symmetry. Furthermore, 
upon applying pressure, the hierarchy of states changes and 
the single stripe state becomes the lowest in energy, as can 
be seen in Fig. [3] thus making fluctuations at q ss = (n, 0) 
the leading mode. 

Discussion 

As mentioned, there are two outstanding experimental 
paradoxes regarding FeSe. The first paradox concerns the 
splitting of the orthorhombic and magnetic transition ob¬ 
served in Fe pnictides, which is taken to an extreme in 
FeSe: the structural transition occurs at T s ~ 90 K, but 
no magnetic order follows. Yet, exactly as in the pnictides, 
DFT calculations reproduce the distorted structure when 
the calculated ground state magnetic structure is used, but 
show no tendency towards orbital ordering or a structural 
distortion if magnetization is kept zero. 

The second paradox deals with the behavior of the critical 
superconducting temperature with pressure T C (P). Typi¬ 
cally, pressure has a tendency to suppress magnetism, so 
in the context of a magnetic pairing mechanism, pressure 
is beneficial to superconductivity when magnetic order is 
present, but it is destructive if it is not. For the nonmag¬ 
netic FeSe, the expectation then is that T c should decrease 
monotonically with pressure. Instead, T c first increases and 
then decreases with pressure, forming a characteristic dome 
shape 37 . In the following we discuss how the J1-J2-J3-K 
model resolves these paradoxes. 

First we analyze the J1-J2-J3-K model parameters given 
in Tablep] and plotted in Fig. [2] The crosses in Figs [2]d , 
[2^, and[2f, show the placement of FeSe at 9 GPa, FeSe at 
0 GPa, and FeTe at 0 GPa, respectively, in the J 1 -J 2 -J 3 - 
K phase diagram. Interestingly, for both FeSe and FeTe 
the calculated ground state at ambient pressure is near a 
phase boundary: between the staggered dimer phase and 
the single stripe phase for FeSe and between the staggered 
dimer phase and double stripe phase for FeTe. Note that 
FeTe appears to be very close to an Ising model because of 
the large K and not because of a large magnetic anisotropy. 

Generally speaking, one can anticipate that, in the ab¬ 
sence of long-range order, spin fluctuations with wave vec¬ 
tors corresponding to the lowest energy states will occur: 


Thus, in FeTe one expects fluctuations with qd s , q<ji, and 
q tr i. None of those would support s± superconductivity 
since only fluctuations with q ~ q ss can pair electrons 
in the standard s± superconducting state. They all break 
tetragonal symmetry, but in different ways, incompatible 
with each other, and cannot all support the same nematic 
state. In FeSe, by contrast, one expects fluctuations with 
q = ( 7 r, Q), where Q = 0, 7r/4, 7 t/ 3, and n/2 (while we 
cannot check this, likely all fluctuations with q = (n,Q), 
where —n/2 < Q < n/2 are supported, cf. Fig. [ 4 ] a). 
This is very different from Fe pnictides, where the single 
stripe state is much lower in energy than all other patterns, 
and therefore the q ss fluctuations dominate. Note that the 
above results rely upon the fact that Fe in FeBS has a large 
local moment (even larger than in DFT),® and cannot be 
obtained by linear response calculations based on a param¬ 
agnetic phas^®. 

The most important consequence of our findings is that 
different spin fluctuations in FeSe (but not FeTe), while 
mutually incompatible with regards to long range magnetic 
order, break tetragonal symmetry in the same way (and the 
same is true for all q = (7 r, Q), — 7t/2 < Q < tt/2). In other 
words, one can have a suppression of long-range magnetic 
ordering due to the competing fluctuations with q = (7r, Q) 
for different Qs, but at the same time these fluctuations all 
break the x y symmetry and do not compete in terms 
of nematicity. Note that the double-stripe fluctuations with 
< 4 ds = (tt/2, tt/2) break a different symmetry, x+y -n- x—y, 
and thus do compete nematically with the single-stripe ones 
in FeTe. Therefore FeSe represents a special case where sev¬ 
eral different types of spin fluctuations are simultaneously 
excited, which prevents them from condensing at any one 
wave vector and forming long range magnetic order, but 
does not prevent the formation of the nematic orthorhom¬ 
bic order. We emphasize that this nematic order, just as 
the underlying incipient magnetic one, is accompanied by 
considerable orbital ordering. We find (see Fig.[4]c) that all 
investigated q = ( 7 r, Q) states induce population imbalance 
between the F e(d xz ) and Fe(d yz ) orbitals on each Fe site of 
the order of (n xz — n yz )/(n xz + n yz ) ~ 8%. This observa¬ 
tion proves that the orbital ordering is not sensitive to the 
magnetic long range order, but only to the nematic order, 
and has multiple ramifications. The orbital ordering can 
be probed experimentally, and was observed in the nematic 
phase at T < T s by the Knight shift anisotropy 3 -, while 
a divergence in 1/TTi, as expected, was only observed at 
much lower temperatures, upon approaching the long range 
magnetic order at T ~ 0. 

Let us now address the intriguing pressure dependence 
of T c . In general, pressure reduces magnetic interactions in 
FeSe. However, the staggered dimer state is suppressed with 
pressure faster than the single stripe state (Fig. [ 3 ]), so that 
instead of multiple competing types of fluctuation we obtain 
a situation similar to the pnictides, where fluctuations with 
q = (tt, 0) decisively dominate. Note that in the s± model 
these are the fluctuations that are responsible for supercon¬ 
ductivity. At ambient pressure, the staggered dimer/trimer 
fluctuations are dominant, but cannot lead to pairing, since 
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the very small FeSe Fermi pockets are not connected by 
q = ( 7 r, Q), where Q ~ n/2. As discussed in Ref. @3] such 
low energy fluctuations with "wrong” momenta are pair¬ 
breaking since they act essentially as impurities (note that 
the situation in FeSe is qualitatively different from previ¬ 
ous discussions in which fluctuations in different channels 
compete®, but can in principle each lead to pairing). 

Under pressure the pairbreaking staggered dimer and 
trimer spin fluctuations are seen to decrease in amplitude 
much more rapidly than the pairing stripe spin fluctuations 
at q = ( 77 ,0). This removal of pairbreaking effects is 
responsible for the initial increase in T c . The further in¬ 
crease of pressure decreases the amplitude of both the pair¬ 
breaking q = ( 7 r, Q) and pairing q = (7r,0) fluctuations, 
leading to the dome-like behavior of T c vs. pressure. 

Even more importantly, the nematic order, which is 
strongest at P = 0, gradually weakens with pressure, as 
the q = (7 t,Q) (Q ~ tt/ 2) fluctuations are suppressed. 
As shown in Fig. [4] b, the density of states at the Fermi 
level N( 0) is strongly decreased in all nematic-compatible 
states compared to the paramagnetic or Neel states, which 
is detrimental for superconductivity (and vice-versa, as ob¬ 
served in BaFe 2 _a;Co. T As 2 !^|). This result indicates that 
long-range orbital, not magnetic, ordering leads to a sharp 
reduction in the Fermi surface and thereby iV(0), which 
is consistent with photoemission and quantum oscillation 
experiment^® 47 . Since T c is exponentially dependent on 
N(0), the suppression of nematicity with pressure is another 
factor ensuring the initial rise of T c . 

Conclusions 

We presented a detailed analysis, based on first principles 
calculations, of magnetic interactions in the FeSe/Te family. 
We show that in FeSe the magnetic interactions are much 
more frustrated than in either FeTe or the Fe pnictides. We 
argue that the simultaneous excitation of spin fluctuations 
with various wave vectors of the type q = (n, Q ) prevent 
long-range magnetic ordering in FeSe, but does allow for 
the usual spin-nematic order accompanied by a ferro-orbital 
order. At zero pressure the leading fluctuations are non¬ 
pairing (in the s± channel) Q = 7 t/2 ones, but pairing 
fluctuations at Q = 0 become the leading fluctuations with 
pressure, which explains the unusual nonmonotonic pressure 
dependence of T c . 

To be able to analyze the emerging situation on a model 
level, we mapped the low-energy energetics onto a three 
neighbors Heisenberg + biquadratic exchange Hamiltonian, 
which we have solved analytically at T = 0 in the mean 
field approximation. It appears that the biquadratic inter¬ 
action is essential to stabilize the observed double stripe 
phase in FeTe; without the extra term, this phase can never 
be the ground state at any choice of parameters. The same 
is true for the staggered dimer phase found to be the DFT 
ground state in FeSe. A nontrivial combination of the bi¬ 
quadratic and third-neighbor exchanges, in addition to the 
usually considered first and second neighbor Heisenberg in¬ 
teractions, ensures the anomalously large splitting of the 
nematic and antiferromagnetic transitions (in FeSe, it leads 
to a total suppression of magnetic ordering). We believe 


that this new perspective on the unusual magnetic physics 
of Fe chalcogenides will be crucial to an explanation of their 
remarkable properties, including perhaps high temperature 
superconductivity in the monolayer FeSe system. 

Methods 

We employed density functional theory and made use 
of three separate, full potential (all electron) codes, ELK, 
WIEN2K, and fplo to calculate the energies. The gener¬ 
alized gradient approximation was used for the exchange- 
correlation functional. We checked for convergence with 
respect to k-points and, for ELK, the number of empty 
states. We calculated the energies of multiple collinear con¬ 
figurations using all three codes for comparison purposes, 
while noncollinear calculations were handled exclusively by 
the ELK code. The comparison of all the different collinear 
configuration energies can be found in the Supplementary 
Materials. 

We used the tetragonal P4/nmm space group (origin 
choice 2) for the crystal structure of FeSe and FeTe in all 
our calculations. The Fe and chalcogenide (Se/Te) ions oc¬ 
cupy the 2a and 2c Wyckoff positions, respectively. The 
lattice parameters for the different materials (and for FeSe, 
under different pressures) are summarized in Table [TJ] We 
note that at low temperatures FeSe is strictly an orthorhom¬ 
bic structure, but this distortion is small and omitting it 
leads to a small magnetoelastic error when compared with 
the exchange parameter energy scales. Furthermore, we are 
interested in the physics that emerges from spin fluctua¬ 
tions that originate in the tetragonal phase. Therefore, for 
FeSe, we defined a volume-conserving effective parameter 
a* = \fah , where a and b are the orthorhombic parameters 
taken from experiment. 

We fit to the Hamiltonian in Equation (JT]) in the usual 
way. The details of how the fit was performed are given in 
the Supplementary Materials. It was not possible to achieve 
a fit that accurately reproduced all energies for all possible 
collinear configurations, so we defined a set of criteria for 
our fitting procedure. The criteria were (1) collinear ground 
states of the J 1 -J 2 -J 3 -K model should be included, (2) low 
energy structures that do not suffer from moment collapse 
under pressure (for FeSe) should be included, (3) local mo¬ 
ments of included structures should be similar, and (4) we 
exclude configurations that yield fits that do not reproduce 
the density functional theory energy hierarchy of the low¬ 
est energy configurations. The fourth criterion is necessary 
because we cannot produce an accurate fit for all configura¬ 
tions, so we decide which features of the density functional 
theory set of energies is important from the point of view 
of fluctuations and frustration, which are the lowest energy 
ones. Given these criteria, we perform the fitting procedure 
using the energies summarized in Fig. [3] 
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SUPPLEMENTARY MATERIAL 


Supplementary Methods 

We calculated the energy of a variety of different 
collinear structures using the three different codes, ELK 4 ®, 
WIEN2k^, and fplg^. The generalized gradient approxi¬ 
mation was used for the exchange-correlation functional. 
The structures in Fig. [5] summarize all of the different con¬ 
figurations that we considered. In Fig. [6] are the energies we 
calculated using these codes. Note that we did not calcu¬ 
late the energy of every configuration using all three codes, 
but there are several points of comparison. For all configura¬ 
tions for which we can make a comparison, there is excellent 
agreement across codes. The most important result of this 
comparison is that there is no ambiguity as to the energy 
hierarchy of the low-lying energy states, it is the same for 
all three codes. We also note that the energy range for 
the different configurations is quite large for both FeSe and 
FeTe, on the scale of 100 — 300 meV for FeSe and 50 — 100 
meV for FeTe. 

We fitted to the Heisenberg model with the ordinary least 
squares (OLS) method using the Heisenberg model coeffi¬ 
cients reported in Table |TTT] for our collinear fits and the 
expression A E{6) = E{8) — E{ 0) = 2K sin 2 (0), see RefP® 
for the definition of 9, for our noncollinear fits. The non- 
collinear energies and the corresponding fits are shown in 

Fig- 13 

As reported in the main article, the collinear fits are very 
good for the included configurations, but there are devia¬ 
tions if we apply the model to configurations excluded due 
to the criteria we outlined in the main Methods section. 
This is a consequence of the itinerant nature of the mag¬ 
netism, which in general cannot be mapped onto a pairwise 
interaction model. We also note that the lower symme¬ 
try magnetic structures, such as those with generic names 
such as "dduuduuu,” suffered moment collapse in FeSe un¬ 
der pressure. The noncollinear fits, on the other hand, are 
excellent. 

It is worth noting that for itinerant magnets the exchange 
model could be potentially improved using an approach sim¬ 
ilar to Moriya 52 and allowing the moment amplitudes to vary 
and by also including Stoner-like onsite terms. We tried in¬ 
cluding terms like this to see how it affected the quality 
of our fits. We found that including these terms does not 
change the fitting results in any qualitative way when using 


the configurations in Table III Furthermore, it did not allow 


us to extend the fit to also reproduce the high-energy con¬ 
figurations from Fig. [5] It is possible, however, that these 
modifications would be important for fluctuations above the 
Neel temperature. 


Phase boundaries of Jj -J r J A -K model 

Here we give additional details of the analytic solution 
of the J 1 -J 2 -J 3 -K model. First, let us ignore the I\ term 
and refresh what is known about the J\ — J2 — J3 Heisen¬ 
berg models. The Ising model has four phases, the checker¬ 
board (cb) phase in Fig. [5^, the double stripe (ds) phase 
in Fig. [5]d, the single stripe (ss) phase in Fig. [5^, and the 
staggered dimers (di) phase in Fig. [bjj The Heisenberg 
model also has four phases, but neither the ds or di phase 
are ground states in the phase diagram. Instead, the four 
phases are the aforementioned cb and ss phases, and in 
addition two spiral phases with spins rotating away from 
the origin as a = n x q x + n y q y , the first with wavevector 
qi = ( 7 r, Q) and the second with q 2 = (Q,Q) (see, e.g., 
RefP®) . Note that at qi(Q — > 0) = (zr , 0), which is the 
ss phase, and at qi(Q —> ir) = q 2 (<3 —t 7 r) = (71-, 7 r), 
which is the cb phase. In both phases Q depends on the 
exchange parameters: Q = cos -1 [( 2 J 2 — J\) /4J3] for qi 
and Q = n— cos -1 [J\/ ( 2 .J 2 + 4 J 3 )] for q 2 . Finally, the an¬ 
alytic expressions for the phase boundaries are summarized 
in Table [TV] 

Adding in the biquadratic term —K (irq • my ) 2 restores 
the ds and di configurations to the phase diagram. The 
allowed wavevectors in the spin spiral phases also become 
dependent on K: Q = cos -1 [(2J 2 — J\) / (4J 3 — 2 1\)\ for 
qi and Q = n — cos -1 [J\/ ( 2 J 2 + 4 J 3 — 2 K)\ for q 2 . The 
analytic expressions for the phase boundaries also change 
and many become A'-dependent as summarized in the last 
column of Table [TV] As K grows so do the areas of stability 
of the ds and di phases. Once K > J\/2, the phase diagram 
becomes indistinguishable from the Ising model. 

FeSeo.sTeo.s 

A notable omission to our results is the case of 
FeSeo. 5 oTeo. 5 o- Unlike the other materials, the structure of 
FeSe 0 . 50 Te 0.50 is not well-defined. A common approach is 
to use lattice parameters from experiment and then choose 
the chalcogenide to be either pure Se or pure Te, assuming 
that the change in the lattice parameters drives the relevant 
physics, such as inducing superconductivity. A check of this 
reveals that this is not entirely the case; there is a signif¬ 
icant energy splitting of the checkerboard, double stripe, 
and zig-zag configurations when Se/Te are swapped, and 
the energy splits are not in the same direction. Using Te 
lowers the checkerboard energy, while it increases the dou¬ 
ble stripe energy, for example. Furthermore, careful exper¬ 
imental analysis reveals that FeSe 0 . 50 Te 0.50 is a disordered 
structure with different heights for Se and Ti^S. Taking this 
into account requires an expensive and non-trivial averag¬ 
ing procedure. While a description of FeSeo. 50 Teo .50 would 
be useful, we put the question aside for now due to the 
complexity of the structure. 
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a single stripe b double stripe C Neel 



q=(7t,0) q=(jt/2,jc/2) q=(jt,jt) 


d staggered dimer e staggered trimer 



q=(jc,jt/2) q=(n,jt/3) 


FIG. 1 . Collinear magnetic structures used for fitting to the J3-J2-J3-K model, (a) single stripe, (b) double stripe, 
(c) checkerboard (Neel), (d) staggered dimer, and (e) staggered trimer state. 


a Ising model b Heisenberg model c K=0A 



J 2 U 1 1 0 05 J 2 U , 1 0 05 J 2 IJ, ' 


FIG. 2 . Classical mean-field phase diagrams, a, The J1-J2-J3 Ising model, b, The J1-J2-J3 Heisenberg model, c, The 
J1-J2-J3-K model with K = 0 . 1 . d, The J1-J2-J3-K model with K = 0.20 where the cross corresponds to FeSe at 9 GPa of 
pressure, e, The J1-J2-J3-K model with K = 0.25 where the cross corresponds to FeSe at 0 GPa of pressure, f, The J1-J2-J3- 
model with K = 0.39 where the cross corresponds to FeTe at 0 GPa of pressure. The length of the cross’ bars in panels 
through f indicate the uncertainty of the fit to the J3-J2-J3-K model. 
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q=( 0 , 7 t) (ji/ 3 ,jt) (71/2,71) (71,71) (71/2, 7 i/ 2 ) 



FIG. 3. Energies of collinear magnetic configurations. The total calculated DFT energies of the collinear magnetic 
configurations used in the fits of the J1-J2-J3-K model. 



FIG. 4. Energies, density of states at ef and orbital order of FeSe at 0 GPa pressure, a The energies of collinear 
magnetic configurations of FeSe at ambient pressure plotted as a function of q. b N(ef) for the lowest energy magnetic 
configurations compared to the nonmagnetic states indicates very small Fermi surfaces in fluctuating magnetic states, c Ferro- 
orbital order in Fe 3 d xz /d yz orbitals measured for several magnetic states. 
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Material 

Ji 

J 2 

Js 

K 




(me 

■V) 



FeSe (0 GPa) 

123.1 ±6.5 

73.0 ±3.3 

18.3 ± 1.8 

30.6 ± 

0.4 

FeSe (4 GPa) 

86.9 ±2.4 

51.9± 1.2 

9.7 ±0.6 

15.7 ± 

0.2 

FeSe (9 GPa) 

51.1 ±0.7 

35.4 ±0.3 

4.9 ±0.2 

10.4 ± 

0.1 

FeTe 

50.7 ± 3.6 

42.8 ± 1.8 

24.4 ± 1.0 

19.7 ± 

0.2 


TABLE I. Heisenberg and biquadratic exchange parameters for FeSe/Te. The parameters for FeSe are reported at 
three different pressures, 0 GPa, 4 GPa, and 9 GPa. FeTe is reported at 0 GPa. The reported uncertainties indicate the 
inaccuracy of the fit to the Ji-J 2 -./ 3 -A' model. 


Material 

a (A) 

c (A) 

ZSe 2Te 

FeSe 

3.76976 

5.52122 

0.2688 

FeTe 

3.81362 

6.25381 

0.2829 

FeSe (4 GPa) 

3.6717 

5.1943 

0.2740 

FeSe (9 GPa) 

3.6049 

5.0304 

0.2839 


TABLE II. Crystal parameters for FeSe and FeTe. The structure parameters and Wyckoff positions for FeSe/Te, with 
FeSe reported at three different pressures. The 4 and 9 GPa parameters are from Ref. EH where we defined an effective 
tetragonal lattice parameter a* = y/ab using the inplane orthorhombic lattice parameters a and b. 
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FIG. 5. Collinear magnetic structures a, Checkerboard (cb). b, Ferromagnetic (fm). c, Parallel Stripes (parastr). d, 
Double Stripe (ds). e, Single Stripe (ss). f, Staggered Dimers (di). g, Zig-zag Stripes (zigzag), h, dduuduuu. i, duuuuuuu. j, 
dduuuuuu. k, udduuuuu. 1, duuuduuu. m, ddduuuuu. n, udduduuu. o, ddduduuu. p, plaquette. q, Staggered Trimers (tri). 
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FIG. 6. Comparison of DFT energies The DFT energies were calculated using elk, wien2k, and fplo. a, FeSe. b, FeTe. 
See Fig. [ 5 ] for the different structures. 


Configuration 

Ji 

J 2 

J 3 

Const. 

Checkerboard 

-2 

2 

2 

1 

Single Stripes 

0 

-2 

2 

1 

Double Stripes 

0 

0 

-2 

1 

Dimers 

-1 

0 

0 

1 

Trimers 

-2/3 

-2/3 

2/3 

1 


TABLE III. The Ji, J2, and J3 coefficients. The coefficients obtained by summing over neighbors of the magnetic structures 
used for fitting to the J1-J2-J3-K model. 






















E(6)-E( 0) (meV/Fe) E(6)-E( 0) (meV/Fe) 


12 



FIG. 7 . Energies of noncollinear structures for FeSe and FeTe as a function of the rotation angle. The dashed 
lines are the model fits, a, FeSe at 0 GPa pressure, b, FeSe at 4 GPa pressure, c, FeSe at 9 GPa pressure, d, FeTe. 


Phase boundary 

Ising 

Heisenberg 

J1-J2-J3-K 

cb/di 

cb/ss 

2J3 + 2 J2 — Ji 


2J3 + 2 J2 — Ji 

ss/di 

2^3 — 2J2 + J\ 


4^/3 — 2J2 + Ji — 2K 

di/ds 

2J3 — Ji 


2J3 — Ji 

cb/q = (tt.Q) 


4J3 + 2 J2 — Ji 

4J3 + 2J2 — Ji — 2 K 

ss/q = (n,Q) 


4J3 — 2 J2 + Ji 

4J 3 - 2 J 2 + Ji - 2/f 

cb/q = (Q,Q) 


4J3 + 2 J2 — Ji 

4J 3 + 2 J 2 - Ji - 2 1< 

q = (n,Q)/q = (Q,Q) 
di/q = (7 v,Q) 
di/q = (Q,Q) 

ds/q = (7r,Q) 

ds/q = {Q,Q) 


2 J3 — J2 

2 J 3 -J 2 -K 

AK (2 J 3 - A') - (2 J 2 - Ji) 2 
(Ji -J2- 4J 3 + 3A') 2 - Ji (Ji + J 2 + A') + 2 Jl 
8J3 (2J 3 - Ji) - (2J 2 - J 2 ) 2 - (2A - Ji) 2 + Jl 
8J3A + 4J 2 A - 47\ 2 - ,7/ 


TABLE IV. Analytical solutions for phase boundaries of the Ising, Heisenberg, and J1-J2-J3-K models. 













